
## @knitr BF_SETUP

rm(list=ls())

# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}


if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


## why not clustered standard errors
## https://stats.stackexchange.com/questions/313618/cox-regression-clustered-standard-errors

## https://cran.r-project.org/web/views/Survival.html
## https://stats.oarc.ucla.edu/r/dae/mixed-effects-cox-regression/
## https://grodri.github.io/survival/frailtyr

library(survival)
## library(coxme) 
library(xtable)
library(ggplot2)
library(english)
library(Publish)
library(showtext)
library(extrafont)
library(mice)
library(mitools)
library(sandwich)
library(lmtest)
library(BFpack)
library(tidyr)
library(gridExtra)
library(latex2exp)

Data <- read.csv('gc-jg-project/kampala/data/data-kampala.csv',
                 stringsAsFactors=FALSE)

Data$log.gdp.cap.lag <- log(Data$gdp.cap.lag)
Data$log.pop.lag <- log(Data$pop.lag)

Data$country_id <- as.numeric(factor(Data$country_name))

m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_polyarchy.lag",
                 "v2x_libdem.lag",     
                 "v2x_partipdem.lag",
                 "v2x_delibdem.lag", 
                 "v2x_egaldem.lag",
                 "log.gdp.cap.lag",
                 "log.pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "country_id")], m = m, maxit = 10, seed = 1510)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()
bf_m1.sidea.results <- list()
bf_m1.sideb.results <- list()
bf_m1.milex.results <- list()
bf_m1.polyarchy.results <- list()
bf_m1.gdp.results <- list()
bf_m1.pop.results <- list()
bf_m1.common.law.results <- list()

BF.m1.sidea <- list()
BF.m1.sideb <- list()
BF.m1.polyarchy <- list()
BF.m1.gdp <- list() 
BF.m1.pop <- list()
BF.m1.common.law <- list()
BF.m1.milex <- list()

BF.m2.sidea <- list()
BF.m2.sideb <- list()
BF.m2.polyarchy <- list()
BF.m2.gdp <- list() 
BF.m2.pop <- list()
BF.m2.common.law <- list()
BF.m2.milex <- list()

BF.m3.sidea <- list()
BF.m3.sideb <- list()
BF.m3.polyarchy <- list()
BF.m3.gdp <- list() 
BF.m3.pop <- list()
BF.m3.common.law <- list()
BF.m3.milex <- list()

BF.m4.sidea <- list()
BF.m4.sideb <- list()
BF.m4.polyarchy <- list()
BF.m4.gdp <- list() 
BF.m4.pop <- list()
BF.m4.common.law <- list()
BF.m4.milex <- list()

BF.m5.sidea <- list()
BF.m5.sideb <- list()
BF.m5.polyarchy <- list()
BF.m5.gdp <- list() 
BF.m5.pop <- list()
BF.m5.common.law <- list()
BF.m5.milex <- list()

BF.m6.sidea <- list()
BF.m6.sideb <- list()
BF.m6.polyarchy <- list()
BF.m6.gdp <- list() 
BF.m6.pop <- list()
BF.m6.common.law <- list()
BF.m6.milex <- list()

m1.bf_objects <- list()
m2.bf_objects <- list()
m3.bf_objects <- list()
m4.bf_objects <- list()
m5.bf_objects <- list()
m6.bf_objects <- list()


for (i in 1:m) {
    complete_data <- complete(imp, i)

    complete_data$.imp <- i
    
    merged_data <- merge(complete_data, 
                         Data[,c("country_id","year","country_name",
                                 "kampala.t0","kampala.t1",
                                 "sidea.1945",
                                 "sideb.1945",
                                 "sidea.force.1945",
                                 "sideb.force.1945",
                                 "sidea.no.force.1945",
                                 "sideb.no.force.1945",
                                 "sidea.force.1991",
                                 "sideb.force.1991",
                                 "sidea.no.force.1991",
                                 "sideb.no.force.1991"
                                 )], 
                        by = c("country_id", "year"),  # Note: quoted variable names
                        all.x = TRUE)  # Keep all rows from complete_data
    
    merged_imputations[[i]] <- merged_data

    current_data <- merged_imputations[[i]]

    m1 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
              sidea.1991
            + sideb.1991
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log.gdp.cap.lag
            + log.pop.lag
            + common.law
            + cluster(country_id)
           ,data=current_data)

     m1.vcov_cluster <- vcovCL(m1, cluster = current_data$country_id)

     m1.model_results[[i]] <- list(
        coefficients = coef(m1),
        variance = m1.vcov_cluster,
        n = nrow(current_data)
       )
     
    m2 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                sidea.1945
            + sideb.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log.gdp.cap.lag
            + log.pop.lag
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m2.vcov_cluster <- vcovCL(m2, cluster = current_data$country_id)

    m2.model_results[[i]] <- list(
        coefficients = coef(m2),
        variance = m2.vcov_cluster,
        n = nrow(current_data)
       )
    
    m3 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1991
                + sideb.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log.gdp.cap.lag
                + log.pop.lag
                + common.law
                + cluster(country_name)
               ,data=current_data)

    m3.vcov_cluster <- vcovCL(m3, cluster = current_data$country_id)

    m3.model_results[[i]] <- list(
        coefficients = coef(m3),
        variance = m3.vcov_cluster,
        n = nrow(current_data)
       )
    
    m4 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1945
            + sideb.force.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log.gdp.cap.lag
            + log.pop.lag
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m4.vcov_cluster <- vcovCL(m4, cluster = current_data$country_id)

    m4.model_results[[i]] <- list(
        coefficients = coef(m4),
        variance = m4.vcov_cluster,
        n = nrow(current_data)
       )
    
    m5 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1991
                + sideb.no.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log.gdp.cap.lag
                + log.pop.lag
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m5.vcov_cluster <- vcovCL(m5, cluster = current_data$country_id)

    m5.model_results[[i]] <- list(
        coefficients = coef(m5),
        variance = m5.vcov_cluster,
        n = nrow(current_data)
       )
    
    m6 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1945
                + sideb.no.force.1945
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log.gdp.cap.lag
                + log.pop.lag
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m6.vcov_cluster <- vcovCL(m6, cluster = current_data$country_id)

    m6.model_results[[i]] <- list(
        coefficients = coef(m6),
        variance = m6.vcov_cluster,
        n = nrow(current_data)
       )

    BF.m1.sidea[[i]] <- BF(m1,hypothesis="sidea.1991<0;sidea.1991=0;sidea.1991>0",complement=TRUE)
    BF.m1.sideb[[i]] <- BF(m1,hypothesis="sideb.1991<0;sideb.1991=0;sideb.1991>0",complement=TRUE)
    BF.m1.polyarchy[[i]] <- BF(m1,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m1.gdp[[i]] <- BF(m1,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m1.pop[[i]] <- BF(m1,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m1.common.law[[i]] <- BF(m1,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m1.milex[[i]] <- BF(m1,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

    BF.m2.sidea[[i]] <- BF(m2,hypothesis="sidea.1945<0;sidea.1945=0;sidea.1945>0",complement=TRUE)
    BF.m2.sideb[[i]] <- BF(m2,hypothesis="sideb.1945<0;sideb.1945=0;sideb.1945>0",complement=TRUE)
    BF.m2.polyarchy[[i]] <- BF(m2,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m2.gdp[[i]] <- BF(m2,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m2.pop[[i]] <- BF(m2,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m2.common.law[[i]] <- BF(m2,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m2.milex[[i]] <- BF(m2,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

    BF.m3.sidea[[i]] <- BF(m3,hypothesis="sidea.force.1991<0;sidea.force.1991=0;sidea.force.1991>0",complement=TRUE)
    BF.m3.sideb[[i]] <- BF(m3,hypothesis="sideb.force.1991<0;sideb.force.1991=0;sideb.force.1991>0",complement=TRUE)
    BF.m3.polyarchy[[i]] <- BF(m3,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m3.gdp[[i]] <- BF(m3,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m3.pop[[i]] <- BF(m3,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m3.common.law[[i]] <- BF(m3,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m3.milex[[i]] <- BF(m3,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

    BF.m4.sidea[[i]] <- BF(m4,hypothesis="sidea.force.1945<0;sidea.force.1945=0;sidea.force.1945>0",complement=TRUE)
    BF.m4.sideb[[i]] <- BF(m4,hypothesis="sideb.force.1945<0;sideb.force.1945=0;sideb.force.1945>0",complement=TRUE)
    BF.m4.polyarchy[[i]] <- BF(m4,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m4.gdp[[i]] <- BF(m4,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m4.pop[[i]] <- BF(m4,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m4.common.law[[i]] <- BF(m4,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m4.milex[[i]] <- BF(m4,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

    BF.m5.sidea[[i]] <- BF(m5,hypothesis="sidea.no.force.1991<0;sidea.no.force.1991=0;sidea.no.force.1991>0",complement=TRUE)
    BF.m5.sideb[[i]] <- BF(m5,hypothesis="sideb.no.force.1991<0;sideb.no.force.1991=0;sideb.no.force.1991>0",complement=TRUE)
    BF.m5.polyarchy[[i]] <- BF(m5,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m5.gdp[[i]] <- BF(m5,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m5.pop[[i]] <- BF(m5,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m5.common.law[[i]] <- BF(m5,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m5.milex[[i]] <- BF(m5,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

    BF.m6.sidea[[i]] <- BF(m6,hypothesis="sidea.no.force.1945<0;sidea.no.force.1945=0;sidea.no.force.1945>0",complement=TRUE)
    BF.m6.sideb[[i]] <- BF(m6,hypothesis="sideb.no.force.1945<0;sideb.no.force.1945=0;sideb.no.force.1945>0",complement=TRUE)
    BF.m6.polyarchy[[i]] <- BF(m6,hypothesis="v2x_libdem.lag<0;v2x_libdem.lag=0;v2x_libdem.lag>0",complement=TRUE)
    BF.m6.gdp[[i]] <- BF(m6,hypothesis="log.gdp.cap.lag<0;log.gdp.cap.lag=0;log.gdp.cap.lag>0",complement=TRUE)
    BF.m6.pop[[i]] <- BF(m6,hypothesis="log.pop.lag<0;log.pop.lag=0;log.pop.lag>0",complement=TRUE)
    BF.m6.common.law[[i]] <- BF(m6,hypothesis="common.law<0;common.law=0;common.law>0",complement=TRUE)
    BF.m6.milex[[i]] <- BF(m6,hypothesis="my.mil.cap.lag<0;my.mil.cap.lag=0;my.mil.cap.lag>0",complement=TRUE)

}


m1.sidea.BFM <- list()
m1.sideb.BFM <- list()
m1.polyarchy.BFM <- list()
m1.gdp.BFM <- list()
m1.pop.BFM <- list()
m1.common_law.BFM <- list()
m1.milex.BFM <- list()

m2.sidea.BFM <- list()
m2.sideb.BFM <- list()
m2.polyarchy.BFM <- list()
m2.gdp.BFM <- list()
m2.pop.BFM <- list()
m2.common_law.BFM <- list()
m2.milex.BFM <- list()

m3.sidea.BFM <- list()
m3.sideb.BFM <- list()
m3.polyarchy.BFM <- list()
m3.gdp.BFM <- list()
m3.pop.BFM <- list()
m3.common_law.BFM <- list()
m3.milex.BFM <- list()

m4.sidea.BFM <- list()
m4.sideb.BFM <- list()
m4.polyarchy.BFM <- list()
m4.gdp.BFM <- list()
m4.pop.BFM <- list()
m4.common_law.BFM <- list()
m4.milex.BFM <- list()

m5.sidea.BFM <- list()
m5.sideb.BFM <- list()
m5.polyarchy.BFM <- list()
m5.gdp.BFM <- list()
m5.pop.BFM <- list()
m5.common_law.BFM <- list()
m5.milex.BFM <- list()

m6.sidea.BFM <- list()
m6.sideb.BFM <- list()
m6.polyarchy.BFM <- list()
m6.gdp.BFM <- list()
m6.pop.BFM <- list()
m6.common_law.BFM <- list()
m6.milex.BFM <- list()


for (i in 1:m) {
# List all BF objects (one per variable)
m1.bf_objects[[i]] <- list(
  sidea = BF.m1.sidea[[i]],
  sideb = BF.m1.sideb[[i]],
  polyarchy = BF.m1.polyarchy[[i]],
  gdp = BF.m1.gdp[[i]],
  pop = BF.m1.pop[[i]],
  common_law = BF.m1.common.law[[i]],
  milex = BF.m1.milex[[i]]
)
    m1.sidea.BFM[[i]] <- BF.m1.sidea[[i]]$BFmatrix_confirmatory
    m1.sideb.BFM[[i]] <- BF.m1.sideb[[i]]$BFmatrix_confirmatory
    m1.polyarchy.BFM[[i]] <- BF.m1.polyarchy[[i]]$BFmatrix_confirmatory
    m1.gdp.BFM[[i]] <- BF.m1.gdp[[i]]$BFmatrix_confirmatory
    m1.pop.BFM[[i]] <- BF.m1.pop[[i]]$BFmatrix_confirmatory
    m1.common_law.BFM[[i]] <- BF.m1.common.law[[i]]$BFmatrix_confirmatory
    m1.milex.BFM[[i]] <- BF.m1.milex[[i]]$BFmatrix_confirmatory

    
m2.bf_objects <- list(
  sidea = BF.m2.sidea,
  sideb = BF.m2.sideb,
  polyarchy = BF.m2.polyarchy,
  gdp = BF.m2.gdp,
  pop = BF.m2.pop,
  common_law = BF.m2.common.law,
  milex = BF.m2.milex
)

    m2.sidea.BFM[[i]] <- BF.m2.sidea[[i]]$BFmatrix_confirmatory
    m2.sideb.BFM[[i]] <- BF.m2.sideb[[i]]$BFmatrix_confirmatory
    m2.polyarchy.BFM[[i]] <- BF.m2.polyarchy[[i]]$BFmatrix_confirmatory
    m2.gdp.BFM[[i]] <- BF.m2.gdp[[i]]$BFmatrix_confirmatory
    m2.pop.BFM[[i]] <- BF.m2.pop[[i]]$BFmatrix_confirmatory
    m2.common_law.BFM[[i]] <- BF.m2.common.law[[i]]$BFmatrix_confirmatory
    m2.milex.BFM[[i]] <- BF.m2.milex[[i]]$BFmatrix_confirmatory

    
m3.bf_objects <- list(
  sidea = BF.m3.sidea,
  sideb = BF.m3.sideb,
  polyarchy = BF.m3.polyarchy,
  gdp = BF.m3.gdp,
  pop = BF.m3.pop,
  common_law = BF.m3.common.law,
  milex = BF.m3.milex
)

    m3.sidea.BFM[[i]] <- BF.m3.sidea[[i]]$BFmatrix_confirmatory
    m3.sideb.BFM[[i]] <- BF.m3.sideb[[i]]$BFmatrix_confirmatory
    m3.polyarchy.BFM[[i]] <- BF.m3.polyarchy[[i]]$BFmatrix_confirmatory
    m3.gdp.BFM[[i]] <- BF.m3.gdp[[i]]$BFmatrix_confirmatory
    m3.pop.BFM[[i]] <- BF.m3.pop[[i]]$BFmatrix_confirmatory
    m3.common_law.BFM[[i]] <- BF.m3.common.law[[i]]$BFmatrix_confirmatory
    m3.milex.BFM[[i]] <- BF.m3.milex[[i]]$BFmatrix_confirmatory

    
m4.bf_objects <- list(
  sidea = BF.m4.sidea,
  sideb = BF.m4.sideb,
  polyarchy = BF.m4.polyarchy,
  gdp = BF.m4.gdp,
  pop = BF.m4.pop,
  common_law = BF.m4.common.law,
  milex = BF.m4.milex
)

    m4.sidea.BFM[[i]] <- BF.m4.sidea[[i]]$BFmatrix_confirmatory
    m4.sideb.BFM[[i]] <- BF.m4.sideb[[i]]$BFmatrix_confirmatory
    m4.polyarchy.BFM[[i]] <- BF.m4.polyarchy[[i]]$BFmatrix_confirmatory
    m4.gdp.BFM[[i]] <- BF.m4.gdp[[i]]$BFmatrix_confirmatory
    m4.pop.BFM[[i]] <- BF.m4.pop[[i]]$BFmatrix_confirmatory
    m4.common_law.BFM[[i]] <- BF.m4.common.law[[i]]$BFmatrix_confirmatory
    m4.milex.BFM[[i]] <- BF.m4.milex[[i]]$BFmatrix_confirmatory

    
m5.bf_objects <- list(
  sidea = BF.m5.sidea,
  sideb = BF.m5.sideb,
  polyarchy = BF.m5.polyarchy,
  gdp = BF.m5.gdp,
  pop = BF.m5.pop,
  common_law = BF.m5.common.law,
  milex = BF.m5.milex
)

    m5.sidea.BFM[[i]] <- BF.m5.sidea[[i]]$BFmatrix_confirmatory
    m5.sideb.BFM[[i]] <- BF.m5.sideb[[i]]$BFmatrix_confirmatory
    m5.polyarchy.BFM[[i]] <- BF.m5.polyarchy[[i]]$BFmatrix_confirmatory
    m5.gdp.BFM[[i]] <- BF.m5.gdp[[i]]$BFmatrix_confirmatory
    m5.pop.BFM[[i]] <- BF.m5.pop[[i]]$BFmatrix_confirmatory
    m5.common_law.BFM[[i]] <- BF.m5.common.law[[i]]$BFmatrix_confirmatory
    m5.milex.BFM[[i]] <- BF.m5.milex[[i]]$BFmatrix_confirmatory

    
m6.bf_objects <- list(
  sidea = BF.m6.sidea,
  sideb = BF.m6.sideb,
  polyarchy = BF.m6.polyarchy,
  gdp = BF.m6.gdp,
  pop = BF.m6.pop,
  common_law = BF.m6.common.law,
  milex = BF.m6.milex
)

    m6.sidea.BFM[[i]] <- BF.m6.sidea[[i]]$BFmatrix_confirmatory
    m6.sideb.BFM[[i]] <- BF.m6.sideb[[i]]$BFmatrix_confirmatory
    m6.polyarchy.BFM[[i]] <- BF.m6.polyarchy[[i]]$BFmatrix_confirmatory
    m6.gdp.BFM[[i]] <- BF.m6.gdp[[i]]$BFmatrix_confirmatory
    m6.pop.BFM[[i]] <- BF.m6.pop[[i]]$BFmatrix_confirmatory
    m6.common_law.BFM[[i]] <- BF.m6.common.law[[i]]$BFmatrix_confirmatory
    m6.milex.BFM[[i]] <- BF.m6.milex[[i]]$BFmatrix_confirmatory

}

mean_bf.m1.sidea <-      exp(apply(simplify2array(lapply(m1.sidea.BFM, log)), 1:2, mean))
mean_bf.m1.sideb <-      exp(apply(simplify2array(lapply(m1.sideb.BFM, log)), 1:2, mean))
mean_bf.m1.polyarchy <-  exp(apply(simplify2array(lapply(m1.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m1.gdp <-        exp(apply(simplify2array(lapply(m1.gdp.BFM, log)), 1:2, mean))
mean_bf.m1.pop <-        exp(apply(simplify2array(lapply(m1.pop.BFM, log)), 1:2, mean))
mean_bf.m1.common_law <- exp(apply(simplify2array(lapply(m1.common_law.BFM, log)), 1:2, mean))
mean_bf.m1.milex <-      exp(apply(simplify2array(lapply(m1.milex.BFM, log)), 1:2, mean))

mean_bf.m2.sidea <-      exp(apply(simplify2array(lapply(m2.sidea.BFM, log)), 1:2, mean))
mean_bf.m2.sideb <-      exp(apply(simplify2array(lapply(m2.sideb.BFM, log)), 1:2, mean))
mean_bf.m2.polyarchy <-  exp(apply(simplify2array(lapply(m2.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m2.gdp <-        exp(apply(simplify2array(lapply(m2.gdp.BFM, log)), 1:2, mean))
mean_bf.m2.pop <-        exp(apply(simplify2array(lapply(m2.pop.BFM, log)), 1:2, mean))
mean_bf.m2.common_law <- exp(apply(simplify2array(lapply(m2.common_law.BFM, log)), 1:2, mean))
mean_bf.m2.milex <-      exp(apply(simplify2array(lapply(m2.milex.BFM, log)), 1:2, mean))

mean_bf.m3.sidea <-      exp(apply(simplify2array(lapply(m3.sidea.BFM, log)), 1:2, mean))
mean_bf.m3.sideb <-      exp(apply(simplify2array(lapply(m3.sideb.BFM, log)), 1:2, mean))
mean_bf.m3.polyarchy <-  exp(apply(simplify2array(lapply(m3.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m3.gdp <-        exp(apply(simplify2array(lapply(m3.gdp.BFM, log)), 1:2, mean))
mean_bf.m3.pop <-        exp(apply(simplify2array(lapply(m3.pop.BFM, log)), 1:2, mean))
mean_bf.m3.common_law <- exp(apply(simplify2array(lapply(m3.common_law.BFM, log)), 1:2, mean))
mean_bf.m3.milex <-      exp(apply(simplify2array(lapply(m3.milex.BFM, log)), 1:2, mean))

mean_bf.m4.sidea <-      exp(apply(simplify2array(lapply(m4.sidea.BFM, log)), 1:2, mean))
mean_bf.m4.sideb <-      exp(apply(simplify2array(lapply(m4.sideb.BFM, log)), 1:2, mean))
mean_bf.m4.polyarchy <-  exp(apply(simplify2array(lapply(m4.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m4.gdp <-        exp(apply(simplify2array(lapply(m4.gdp.BFM, log)), 1:2, mean))
mean_bf.m4.pop <-        exp(apply(simplify2array(lapply(m4.pop.BFM, log)), 1:2, mean))
mean_bf.m4.common_law <- exp(apply(simplify2array(lapply(m4.common_law.BFM, log)), 1:2, mean))
mean_bf.m4.milex <-      exp(apply(simplify2array(lapply(m4.milex.BFM, log)), 1:2, mean))

mean_bf.m5.sidea <-      exp(apply(simplify2array(lapply(m5.sidea.BFM, log)), 1:2, mean))
mean_bf.m5.sideb <-      exp(apply(simplify2array(lapply(m5.sideb.BFM, log)), 1:2, mean))
mean_bf.m5.polyarchy <-  exp(apply(simplify2array(lapply(m5.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m5.gdp <-        exp(apply(simplify2array(lapply(m5.gdp.BFM, log)), 1:2, mean))
mean_bf.m5.pop <-        exp(apply(simplify2array(lapply(m5.pop.BFM, log)), 1:2, mean))
mean_bf.m5.common_law <- exp(apply(simplify2array(lapply(m5.common_law.BFM, log)), 1:2, mean))
mean_bf.m5.milex <-      exp(apply(simplify2array(lapply(m5.milex.BFM, log)), 1:2, mean))

mean_bf.m6.sidea <-      exp(apply(simplify2array(lapply(m6.sidea.BFM, log)), 1:2, mean))
mean_bf.m6.sideb <-      exp(apply(simplify2array(lapply(m6.sideb.BFM, log)), 1:2, mean))
mean_bf.m6.polyarchy <-  exp(apply(simplify2array(lapply(m6.polyarchy.BFM, log)), 1:2, mean))
mean_bf.m6.gdp <-        exp(apply(simplify2array(lapply(m6.gdp.BFM, log)), 1:2, mean))
mean_bf.m6.pop <-        exp(apply(simplify2array(lapply(m6.pop.BFM, log)), 1:2, mean))
mean_bf.m6.common_law <- exp(apply(simplify2array(lapply(m6.common_law.BFM, log)), 1:2, mean))
mean_bf.m6.milex <-      exp(apply(simplify2array(lapply(m6.milex.BFM, log)), 1:2, mean))


my.extract_pairwise_bf <- function(bf_obj) {
  m <- bf_obj
  data.frame(
    H1_vs_H2 = m[1,2],  # BF(H1 vs H2)
    H1_vs_H3 = m[1,3],  # BF(H1 vs H3)
    H2_vs_H1 = m[2,1],   # BF(H2 vs H1)
    H2_vs_H3 = m[2,3],  # BF(H2 vs H3)
    H3_vs_H1 = m[3,1],  # BF(H3 vs H1)
    H3_vs_H2 = m[3,2]   # BF(H3 vs H2)
  )
}


m1.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m1.sidea),
    my.extract_pairwise_bf(mean_bf.m1.sideb),
    my.extract_pairwise_bf(mean_bf.m1.polyarchy),
    my.extract_pairwise_bf(mean_bf.m1.gdp),
    my.extract_pairwise_bf(mean_bf.m1.pop),
    my.extract_pairwise_bf(mean_bf.m1.common_law),
    my.extract_pairwise_bf(mean_bf.m1.milex)
)
m1.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m1.bf_df <- m1.pairwise_bfs

m2.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m2.sidea),
    my.extract_pairwise_bf(mean_bf.m2.sideb),
    my.extract_pairwise_bf(mean_bf.m2.polyarchy),
    my.extract_pairwise_bf(mean_bf.m2.gdp),
    my.extract_pairwise_bf(mean_bf.m2.pop),
    my.extract_pairwise_bf(mean_bf.m2.common_law),
    my.extract_pairwise_bf(mean_bf.m2.milex)
)
m2.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m2.bf_df <- m2.pairwise_bfs

m3.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m3.sidea),
    my.extract_pairwise_bf(mean_bf.m3.sideb),
    my.extract_pairwise_bf(mean_bf.m3.polyarchy),
    my.extract_pairwise_bf(mean_bf.m3.gdp),
    my.extract_pairwise_bf(mean_bf.m3.pop),
    my.extract_pairwise_bf(mean_bf.m3.common_law),
    my.extract_pairwise_bf(mean_bf.m3.milex)
)
m3.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m3.bf_df <- m3.pairwise_bfs

m4.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m4.sidea),
    my.extract_pairwise_bf(mean_bf.m4.sideb),
    my.extract_pairwise_bf(mean_bf.m4.polyarchy),
    my.extract_pairwise_bf(mean_bf.m4.gdp),
    my.extract_pairwise_bf(mean_bf.m4.pop),
    my.extract_pairwise_bf(mean_bf.m4.common_law),
    my.extract_pairwise_bf(mean_bf.m4.milex)
)
m4.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m4.bf_df <- m4.pairwise_bfs


m5.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m5.sidea),
    my.extract_pairwise_bf(mean_bf.m5.sideb),
    my.extract_pairwise_bf(mean_bf.m5.polyarchy),
    my.extract_pairwise_bf(mean_bf.m5.gdp),
    my.extract_pairwise_bf(mean_bf.m5.pop),
    my.extract_pairwise_bf(mean_bf.m5.common_law),
    my.extract_pairwise_bf(mean_bf.m5.milex)
)
m5.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m5.bf_df <- m5.pairwise_bfs

m6.pairwise_bfs <- rbind(
    my.extract_pairwise_bf(mean_bf.m6.sidea),
    my.extract_pairwise_bf(mean_bf.m6.sideb),
    my.extract_pairwise_bf(mean_bf.m6.polyarchy),
    my.extract_pairwise_bf(mean_bf.m6.gdp),
    my.extract_pairwise_bf(mean_bf.m6.pop),
    my.extract_pairwise_bf(mean_bf.m6.common_law),
    my.extract_pairwise_bf(mean_bf.m6.milex)
)
m6.pairwise_bfs$Variable <- c("sidea",
                              "sideb",
                              "polyarchy",
                              "gdp",
                              "pop",
                              "common_law",
                              "milex")
m6.bf_df <- m6.pairwise_bfs




m1.bf_long <- pivot_longer(
  m1.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m1.bf_long <- pivot_longer(
  m1.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m2.bf_long <- pivot_longer(
  m2.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m3.bf_long <- pivot_longer(
  m3.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m4.bf_long <- pivot_longer(
  m4.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m5.bf_long <- pivot_longer(
  m5.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)

m6.bf_long <- pivot_longer(
  m6.bf_df,
  cols = c(H1_vs_H2, H1_vs_H3, H2_vs_H1, H2_vs_H3,
           H3_vs_H1, H3_vs_H2
           ),
  names_to = "Comparison",
  values_to = "BF"
)


# Log-transform for better visualization
m1.bf_long$log_BF <- log10(m1.bf_long$BF)
m2.bf_long$log_BF <- log10(m2.bf_long$BF)
m3.bf_long$log_BF <- log10(m3.bf_long$BF)
m4.bf_long$log_BF <- log10(m4.bf_long$BF)
m5.bf_long$log_BF <- log10(m5.bf_long$BF)
m6.bf_long$log_BF <- log10(m6.bf_long$BF)



variable_order <- c(
  "sidea",
  "sideb",
  "milex",
  "polyarchy",
  "gdp",
  "pop",
  "common_law"
  )

m1.bf_long$Variable.fct <- factor(
  m1.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m2.bf_long$Variable.fct <- factor(
  m2.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m3.bf_long$Variable.fct <- factor(
  m3.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m4.bf_long$Variable.fct <- factor(
  m1.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m5.bf_long$Variable.fct <- factor(
  m5.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m6.bf_long$Variable.fct <- factor(
  m6.bf_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)




custom_labels <- c(
  "H1_vs_H2" = "$H_1: \\beta$ < 0 vs. $H_0: \\beta$ = 0",   # H1: β < 0 vs. H2: β = 0
  "H1_vs_H3" = "$H_1: \\beta$ < 0 vs. $H_0: \\beta$ > 0",   # H1: β < 0 vs. H3: β > 0
  "H2_vs_H3" = "$H_1: \\beta$ = 0 vs. $H_0: \\beta$ > 0",    # H2: β = 0 vs. H3: β > 0
  "H2_vs_H1" = "$H_1: \\beta$ = 0 vs. $H_0: \\beta$ < 0",   # H1: β < 0 vs. H2: β = 0
  "H3_vs_H1" = "$H_1: \\beta$ > 0 vs. $H_0: \\beta$ < 0",   # H1: β < 0 vs. H3: β > 0
  "H3_vs_H2" = "$H_1: \\beta$ > 0 vs. $H_0: \\beta$ = 0"    # H2: β = 0 vs. H3: β > 0
)



custom_y_labels <- c(
  "sidea" = "MID Side A",
  "sideb" = "MID Side B",
  "polyarchy" = "Polyarchy",
  "gdp" = "GDP cap (log)",
  "pop" = "Population (log)",
  "common_law" = "Common Law",
  "milex" = "Milex cap"
)


n_blocks <- 7  # Number of blocks
base <- c(1:2,4)
sequence <- rep(base, n_blocks) + 6 * rep(0:(n_blocks - 1), each = length(base))
sequence <- sort(sequence)

global_min <- min(c(m1.bf_long$log_BF, m2.bf_long$log_BF,
                    m3.bf_long$log_BF, m4.bf_long$log_BF,
                    m5.bf_long$log_BF, m6.bf_long$log_BF), na.rm = TRUE)
global_max <- max(c(m1.bf_long$log_BF, m2.bf_long$log_BF,
                    m3.bf_long$log_BF, m4.bf_long$log_BF,
                    m5.bf_long$log_BF, m6.bf_long$log_BF), na.rm = TRUE)


# Define shared fill scale
shared_fill_scale <- scale_fill_gradient2(
  low = "blue",
  mid = "white",
  high = "red",
  midpoint = 0,
  name = "log10(BF)",
  limits = c(global_min, global_max)
)


m1.plot <- ggplot(m1.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "",##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Overall, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m2.plot <- ggplot(m2.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
    geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Overall, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m3.plot <- ggplot(m3.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Force, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m4.plot <- ggplot(m4.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Force, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m5.plot <- ggplot(m5.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "No force, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m6.plot <- ggplot(m6.bf_long[sequence,], aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "No force, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )

## knitr BF_FIGURE
bf.figure.notused <- grid.arrange(m1.plot, m3.plot, m5.plot,
             m2.plot, m4.plot, m6.plot, 
             nrow = 2)


bf.1.sidea.1991 <- round(m1.bf_long[1,4],1)
bf.2.sidea.1991 <- round(m1.bf_long[2,4],1)
bf.3.sidea.1991 <- round(m1.bf_long[4,4],1)

bf.1.sideb.1991 <- round(m1.bf_long[7,4],1)
bf.2.sideb.1991 <- round(m1.bf_long[8,4],1)
bf.3.sideb.1991 <- round(m1.bf_long[10,4],1)




## two-tailed hypothesis

post_odds.m1.sidea <- list()
post_odds.m1.sideb <- list()     
post_odds.m1.milex <- list()
post_odds.m1.polyarchy <- list()
post_odds.m1.gdp <- list()       
post_odds.m1.pop <- list()     
post_odds.m1.common.law <- list()

BF.2tail.m1.sidea <- list()
BF.2tail.m1.sideb        <- list()
BF.2tail.m1.milex        <- list()
BF.2tail.m1.polyarchy    <- list()
BF.2tail.m1.gdp          <- list()
BF.2tail.m1.pop          <- list()
BF.2tail.m1.common.law   <- list()

post_odds.m2.sidea       <- list()
post_odds.m2.sideb       <- list()
post_odds.m2.milex       <- list()
post_odds.m2.polyarchy   <- list()
post_odds.m2.gdp         <- list()
post_odds.m2.pop         <- list()
post_odds.m2.common.law  <- list()

BF.2tail.m2.sidea        <- list()
BF.2tail.m2.sideb        <- list()
BF.2tail.m2.milex        <- list()
BF.2tail.m2.polyarchy    <- list()
BF.2tail.m2.gdp          <- list()
BF.2tail.m2.pop          <- list()
BF.2tail.m2.common.law   <- list()

post_odds.m3.sidea       <- list()
post_odds.m3.sideb       <- list()
post_odds.m3.milex       <- list()
post_odds.m3.polyarchy   <- list()
post_odds.m3.gdp         <- list()
post_odds.m3.pop         <- list()
post_odds.m3.common.law  <- list()

BF.2tail.m3.sidea        <- list()
BF.2tail.m3.sideb        <- list()
BF.2tail.m3.milex        <- list()
BF.2tail.m3.polyarchy    <- list()
BF.2tail.m3.gdp          <- list()
BF.2tail.m3.pop          <- list()
BF.2tail.m3.common.law   <- list()

post_odds.m4.sidea       <- list()
post_odds.m4.sideb       <- list()
post_odds.m4.milex       <- list()
post_odds.m4.polyarchy   <- list()
post_odds.m4.gdp         <- list()
post_odds.m4.pop         <- list()
post_odds.m4.common.law  <- list()

BF.2tail.m4.sidea        <- list()
BF.2tail.m4.sideb        <- list()
BF.2tail.m4.milex        <- list()
BF.2tail.m4.polyarchy    <- list()
BF.2tail.m4.gdp          <- list()
BF.2tail.m4.pop          <- list()
BF.2tail.m4.common.law   <- list()

post_odds.m5.sidea       <- list()
post_odds.m5.sideb       <- list()
post_odds.m5.milex       <- list()
post_odds.m5.polyarchy   <- list()
post_odds.m5.gdp         <- list()
post_odds.m5.pop         <- list()
post_odds.m5.common.law  <- list()

BF.2tail.m5.sidea        <- list()
BF.2tail.m5.sideb        <- list()
BF.2tail.m5.milex        <- list()
BF.2tail.m5.polyarchy    <- list()
BF.2tail.m5.gdp          <- list()
BF.2tail.m5.pop          <- list()
BF.2tail.m5.common.law   <- list()

post_odds.m6.sidea       <- list()
post_odds.m6.sideb       <- list()
post_odds.m6.milex       <- list()
post_odds.m6.polyarchy   <- list()
post_odds.m6.gdp         <- list()
post_odds.m6.pop         <- list()
post_odds.m6.common.law  <- list()

BF.2tail.m6.sidea        <- list()
BF.2tail.m6.sideb        <- list()
BF.2tail.m6.milex        <- list()
BF.2tail.m6.polyarchy    <- list()
BF.2tail.m6.gdp          <- list()
BF.2tail.m6.pop          <- list()
BF.2tail.m6.common.law   <- list()

prior_odds <- 2

for (i in 1:m) {

    ## To avoid infinte BF for gdp:
    ## Pr(GDP=0)=0

    if (BF.m1.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m1.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m1.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m2.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m2.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m2.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m3.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m3.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m3.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m4.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m4.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m4.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m5.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m5.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m5.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m6.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m6.gdp[[i]]$PHP_confirmatory[3] == 1) {
        BF.m6.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}


    if (BF.m1.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m1.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m1.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m2.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m2.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m2.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m3.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m3.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m3.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m4.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m4.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m4.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m5.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m5.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m5.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}

    if (BF.m6.gdp[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m6.gdp[[i]]$PHP_confirmatory[3] == .999) {
        BF.m6.gdp[[i]]$PHP_confirmatory[2:3] <- c(0.0001, 0.9999)}


  
    ## to avoid infinite BF for milex
    ## pr(milex=0)=0
    
    if (BF.m1.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m1.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m1.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m2.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m2.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m2.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m3.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m3.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m3.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m4.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m4.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m4.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m5.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m5.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m5.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m6.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m6.milex[[i]]$PHP_confirmatory[1] == 1) {
        BF.m6.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    
    if (BF.m1.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m1.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m1.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m2.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m2.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m2.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m3.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m3.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m3.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m4.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m4.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m4.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m5.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m5.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m5.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}

    if (BF.m6.milex[[i]]$PHP_confirmatory[2] == 0 && 
        BF.m6.milex[[i]]$PHP_confirmatory[1] == .999) {
        BF.m6.milex[[i]]$PHP_confirmatory[1:2] <- c(0.9999,0.0001)}


    
    
post_odds.m1.sidea[[i]]      <- sum(BF.m1.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.sidea[[i]]$PHP_confirmatory[2]
post_odds.m1.sideb[[i]]      <- sum(BF.m1.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.sideb[[i]]$PHP_confirmatory[2]
post_odds.m1.milex[[i]]      <- sum(BF.m1.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.milex[[i]]$PHP_confirmatory[2]
post_odds.m1.polyarchy[[i]]  <- sum(BF.m1.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m1.gdp[[i]]        <- sum(BF.m1.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.gdp[[i]]$PHP_confirmatory[2]
post_odds.m1.pop[[i]]        <- sum(BF.m1.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.pop[[i]]$PHP_confirmatory[2]
post_odds.m1.common.law[[i]] <- sum(BF.m1.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m1.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m1.sidea[[i]]      <- post_odds.m1.sidea[[i]]/prior_odds
BF.2tail.m1.sideb[[i]]      <- post_odds.m1.sideb[[i]]/prior_odds    
BF.2tail.m1.milex[[i]]      <- post_odds.m1.milex[[i]]/prior_odds    
BF.2tail.m1.polyarchy[[i]]  <- post_odds.m1.polyarchy[[i]]/prior_odds
BF.2tail.m1.gdp[[i]]        <- post_odds.m1.gdp[[i]]/prior_odds
BF.2tail.m1.pop[[i]]        <- post_odds.m1.pop[[i]]/prior_odds      
BF.2tail.m1.common.law[[i]] <- post_odds.m1.common.law[[i]]/prior_odds

post_odds.m2.sidea[[i]] <- sum(BF.m2.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.sidea[[i]]$PHP_confirmatory[2]
post_odds.m2.sideb[[i]] <- sum(BF.m2.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.sideb[[i]]$PHP_confirmatory[2]
post_odds.m2.milex[[i]] <- sum(BF.m2.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.milex[[i]]$PHP_confirmatory[2]
post_odds.m2.polyarchy[[i]] <- sum(BF.m2.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m2.gdp[[i]] <- sum(BF.m2.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.gdp[[i]]$PHP_confirmatory[2]
post_odds.m2.pop[[i]] <- sum(BF.m2.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.pop[[i]]$PHP_confirmatory[2]
post_odds.m2.common.law[[i]] <- sum(BF.m2.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m2.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m2.sidea[[i]] <- post_odds.m2.sidea[[i]]/prior_odds
BF.2tail.m2.sideb[[i]] <- post_odds.m2.sideb[[i]]/prior_odds    
BF.2tail.m2.milex[[i]] <- post_odds.m2.milex[[i]]/prior_odds    
BF.2tail.m2.polyarchy[[i]] <- post_odds.m2.polyarchy[[i]]/prior_odds
BF.2tail.m2.gdp[[i]] <- post_odds.m2.gdp[[i]]/prior_odds
BF.2tail.m2.pop[[i]] <- post_odds.m2.pop[[i]]/prior_odds      
BF.2tail.m2.common.law[[i]] <- post_odds.m2.common.law[[i]]/prior_odds


post_odds.m3.sidea[[i]] <- sum(BF.m3.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.sidea[[i]]$PHP_confirmatory[2]
post_odds.m3.sideb[[i]] <- sum(BF.m3.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.sideb[[i]]$PHP_confirmatory[2]
post_odds.m3.milex[[i]] <- sum(BF.m3.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.milex[[i]]$PHP_confirmatory[2]
post_odds.m3.polyarchy[[i]] <- sum(BF.m3.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m3.gdp[[i]] <- sum(BF.m3.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.gdp[[i]]$PHP_confirmatory[2]
post_odds.m3.pop[[i]] <- sum(BF.m3.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.pop[[i]]$PHP_confirmatory[2]
post_odds.m3.common.law[[i]] <- sum(BF.m3.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m3.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m3.sidea[[i]] <- post_odds.m3.sidea[[i]]/prior_odds
BF.2tail.m3.sideb[[i]] <- post_odds.m3.sideb[[i]]/prior_odds    
BF.2tail.m3.milex[[i]] <- post_odds.m3.milex[[i]]/prior_odds    
BF.2tail.m3.polyarchy[[i]] <- post_odds.m3.polyarchy[[i]]/prior_odds
BF.2tail.m3.gdp[[i]] <- post_odds.m3.gdp[[i]]/prior_odds
BF.2tail.m3.pop[[i]] <- post_odds.m3.pop[[i]]/prior_odds      
BF.2tail.m3.common.law[[i]] <- post_odds.m3.common.law[[i]]/prior_odds


post_odds.m4.sidea[[i]] <- sum(BF.m4.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.sidea[[i]]$PHP_confirmatory[2]
post_odds.m4.sideb[[i]] <- sum(BF.m4.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.sideb[[i]]$PHP_confirmatory[2]
post_odds.m4.milex[[i]] <- sum(BF.m4.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.milex[[i]]$PHP_confirmatory[2]
post_odds.m4.polyarchy[[i]] <- sum(BF.m4.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m4.gdp[[i]] <- sum(BF.m4.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.gdp[[i]]$PHP_confirmatory[2]
post_odds.m4.pop[[i]] <- sum(BF.m4.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.pop[[i]]$PHP_confirmatory[2]
post_odds.m4.common.law[[i]] <- sum(BF.m4.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m4.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m4.sidea[[i]] <- post_odds.m4.sidea[[i]]/prior_odds
BF.2tail.m4.sideb[[i]] <- post_odds.m4.sideb[[i]]/prior_odds    
BF.2tail.m4.milex[[i]] <- post_odds.m4.milex[[i]]/prior_odds    
BF.2tail.m4.polyarchy[[i]] <- post_odds.m4.polyarchy[[i]]/prior_odds
BF.2tail.m4.gdp[[i]] <- post_odds.m4.gdp[[i]]/prior_odds
BF.2tail.m4.pop[[i]] <- post_odds.m4.pop[[i]]/prior_odds      
BF.2tail.m4.common.law[[i]] <- post_odds.m4.common.law[[i]]/prior_odds


post_odds.m5.sidea[[i]] <- sum(BF.m5.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.sidea[[i]]$PHP_confirmatory[2]
post_odds.m5.sideb[[i]] <- sum(BF.m5.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.sideb[[i]]$PHP_confirmatory[2]
post_odds.m5.milex[[i]] <- sum(BF.m5.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.milex[[i]]$PHP_confirmatory[2]
post_odds.m5.polyarchy[[i]] <- sum(BF.m5.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m5.gdp[[i]] <- sum(BF.m5.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.gdp[[i]]$PHP_confirmatory[2]
post_odds.m5.pop[[i]] <- sum(BF.m5.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.pop[[i]]$PHP_confirmatory[2]
post_odds.m5.common.law[[i]] <- sum(BF.m5.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m5.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m5.sidea[[i]] <- post_odds.m5.sidea[[i]]/prior_odds
BF.2tail.m5.sideb[[i]] <- post_odds.m5.sideb[[i]]/prior_odds    
BF.2tail.m5.milex[[i]] <- post_odds.m5.milex[[i]]/prior_odds    
BF.2tail.m5.polyarchy[[i]] <- post_odds.m5.polyarchy[[i]]/prior_odds
BF.2tail.m5.gdp[[i]] <- post_odds.m5.gdp[[i]]/prior_odds
BF.2tail.m5.pop[[i]] <- post_odds.m5.pop[[i]]/prior_odds      
BF.2tail.m5.common.law[[i]] <- post_odds.m5.common.law[[i]]/prior_odds


post_odds.m6.sidea[[i]] <- sum(BF.m6.sidea[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.sidea[[i]]$PHP_confirmatory[2]
post_odds.m6.sideb[[i]] <- sum(BF.m6.sideb[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.sideb[[i]]$PHP_confirmatory[2]
post_odds.m6.milex[[i]] <- sum(BF.m6.milex[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.milex[[i]]$PHP_confirmatory[2]
post_odds.m6.polyarchy[[i]] <- sum(BF.m6.polyarchy[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.polyarchy[[i]]$PHP_confirmatory[2]
post_odds.m6.gdp[[i]] <- sum(BF.m6.gdp[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.gdp[[i]]$PHP_confirmatory[2]
post_odds.m6.pop[[i]] <- sum(BF.m6.pop[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.pop[[i]]$PHP_confirmatory[2]
post_odds.m6.common.law[[i]] <- sum(BF.m6.common.law[[i]]$PHP_confirmatory[c(1,3)])/BF.m6.common.law[[i]]$PHP_confirmatory[2]

BF.2tail.m6.sidea[[i]] <- post_odds.m6.sidea[[i]]/prior_odds
BF.2tail.m6.sideb[[i]] <- post_odds.m6.sideb[[i]]/prior_odds    
BF.2tail.m6.milex[[i]] <- post_odds.m6.milex[[i]]/prior_odds    
BF.2tail.m6.polyarchy[[i]] <- post_odds.m6.polyarchy[[i]]/prior_odds
BF.2tail.m6.gdp[[i]] <- post_odds.m6.gdp[[i]]/prior_odds
BF.2tail.m6.pop[[i]] <- post_odds.m6.pop[[i]]/prior_odds      
BF.2tail.m6.common.law[[i]] <- post_odds.m6.common.law[[i]]/prior_odds
}


bf.2tail_long <- rbind(
    BF.2tail.m1.sidea,
    BF.2tail.m2.sidea,
    BF.2tail.m3.sidea,
    BF.2tail.m4.sidea,
    BF.2tail.m5.sidea,
    BF.2tail.m6.sidea,
    BF.2tail.m1.sideb,
    BF.2tail.m2.sideb,
    BF.2tail.m3.sideb,
    BF.2tail.m4.sideb,
    BF.2tail.m5.sideb,
    BF.2tail.m6.sideb,
    BF.2tail.m1.milex,
    BF.2tail.m2.milex,
    BF.2tail.m3.milex,
    BF.2tail.m4.milex,
    BF.2tail.m5.milex,
    BF.2tail.m6.milex,
    BF.2tail.m1.polyarchy, 
    BF.2tail.m2.polyarchy, 
    BF.2tail.m3.polyarchy, 
    BF.2tail.m4.polyarchy, 
    BF.2tail.m5.polyarchy, 
    BF.2tail.m6.polyarchy, 
    BF.2tail.m1.gdp,
    BF.2tail.m2.gdp,
    BF.2tail.m3.gdp,
    BF.2tail.m4.gdp,
    BF.2tail.m5.gdp,
    BF.2tail.m6.gdp,
    BF.2tail.m1.pop,
    BF.2tail.m2.pop,
    BF.2tail.m3.pop,
    BF.2tail.m4.pop,
    BF.2tail.m5.pop,
    BF.2tail.m6.pop,
    BF.2tail.m1.common.law,
    BF.2tail.m2.common.law,
    BF.2tail.m3.common.law,
    BF.2tail.m4.common.law,
    BF.2tail.m5.common.law,
    BF.2tail.m6.common.law)


bf_matrix.2tail <- as.matrix(bf.2tail_long)
bf_matrix.2tail <- matrix(as.numeric(bf_matrix.2tail), nrow = nrow(bf_matrix.2tail))
rownames(bf_matrix.2tail) <- rownames(bf.2tail_long)  # Preserve row names

# Compute geometric mean by row
geo_mean_rows <- apply(bf_matrix.2tail, 1, function(x) {
  exp(mean(log(x)))
})

# Create results table
bf.2tail_long1 <- data.frame(
  RowName = rownames(bf_matrix.2tail),
  GeometricMean = geo_mean_rows
)

bf.2tail_long1$log_BF <- log10(bf.2tail_long1$GeometricMean)

m1.bf.2tail_long <- bf.2tail_long1[seq(1,42,6),3]
m2.bf.2tail_long <- bf.2tail_long1[seq(2,42,6),3]
m3.bf.2tail_long <- bf.2tail_long1[seq(3,42,6),3]
m4.bf.2tail_long <- bf.2tail_long1[seq(4,42,6),3]
m5.bf.2tail_long <- bf.2tail_long1[seq(5,42,6),3]
m6.bf.2tail_long <- bf.2tail_long1[seq(6,42,6),3]


m1.bf.2tail_long <- data.frame('log_BF'=m1.bf.2tail_long)

m1.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m1.bf.2tail_long$Variable.fct <- factor(
  m1.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m1.bf.2tail_long$BF <- exp(m1.bf.2tail_long$log_BF)

m1.bf.2tail_long$Comparison <- "H4_vs_H2"


## [c(1:2,15,
##    3:4,16,
##    5:6,17,
##    7:8,18,
##    9:10,19,
##    11:12,20,
##    13:14,21),]


m2.bf.2tail_long <- data.frame('log_BF'=m2.bf.2tail_long)

m2.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m2.bf.2tail_long$Variable.fct <- factor(
  m2.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m2.bf.2tail_long$BF <- exp(m2.bf.2tail_long$log_BF)

m2.bf.2tail_long$Comparison <- "H4_vs_H2"


m3.bf.2tail_long <- data.frame('log_BF'=m3.bf.2tail_long)

m3.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m3.bf.2tail_long$Variable.fct <- factor(
  m3.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m3.bf.2tail_long$BF <- exp(m3.bf.2tail_long$log_BF)

m3.bf.2tail_long$Comparison <- "H4_vs_H2"


m4.bf.2tail_long <- data.frame('log_BF'=m4.bf.2tail_long)

m4.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m4.bf.2tail_long$Variable.fct <- factor(
  m4.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m4.bf.2tail_long$BF <- exp(m4.bf.2tail_long$log_BF)

m4.bf.2tail_long$Comparison <- "H4_vs_H2"


m5.bf.2tail_long <- data.frame('log_BF'=m5.bf.2tail_long)

m5.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m5.bf.2tail_long$Variable.fct <- factor(
  m5.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m5.bf.2tail_long$BF <- exp(m5.bf.2tail_long$log_BF)

m5.bf.2tail_long$Comparison <- "H4_vs_H2"


m6.bf.2tail_long <- data.frame('log_BF'=m6.bf.2tail_long)

m6.bf.2tail_long$Variable <- c("sidea",
                               "sideb",
                               "milex",
                               "polyarchy",
                               "gdp",
                               "pop",
                               "common_law")

m6.bf.2tail_long$Variable.fct <- factor(
  m6.bf.2tail_long$Variable,
  levels = rev(variable_order)  # Assign the predefined order
)

m6.bf.2tail_long$BF <- exp(m6.bf.2tail_long$log_BF)

m6.bf.2tail_long$Comparison <- "H4_vs_H2"


n_blocks <- 7  # Number of blocks
base <- c(1,6)
sequence <- rep(base, n_blocks) + 6 * rep(0:(n_blocks - 1), each = length(base))
sequence <- sort(sequence)

m1.bf.new_long <- rbind(m1.bf_long[sequence,],m1.bf.2tail_long)
m2.bf.new_long <- rbind(m2.bf_long[sequence,],m2.bf.2tail_long)
m3.bf.new_long <- rbind(m3.bf_long[sequence,],m3.bf.2tail_long)
m4.bf.new_long <- rbind(m4.bf_long[sequence,],m4.bf.2tail_long)
m5.bf.new_long <- rbind(m5.bf_long[sequence,],m5.bf.2tail_long)
m6.bf.new_long <- rbind(m6.bf_long[sequence,],m6.bf.2tail_long)


custom_labels2 <- c(
  "H1_vs_H2" = "$H_1: \\beta$ < 0 vs. $H_0: \\beta$ = 0",   # H1: β < 0 vs. H2: β = 0
  "H1_vs_H3" = "$H_1: \\beta$ < 0 vs. $H_0: \\beta$ > 0",   # H1: β < 0 vs. H3: β > 0
  "H2_vs_H3" = "$H_1: \\beta$ = 0 vs. $H_0: \\beta$ > 0",    # H2: β = 0 vs. H3: β > 0
  "H2_vs_H1" = "$H_1: \\beta$ = 0 vs. $H_0: \\beta$ < 0",   # H1: β < 0 vs. H2: β = 0
  "H3_vs_H1" = "$H_1: \\beta$ > 0 vs. $H_0: \\beta$ < 0",   # H1: β < 0 vs. H3: β > 0
  "H3_vs_H2" = "$H_1: \\beta$ > 0 vs. $H_0: \\beta$ = 0",    # H2: β = 0 vs. H3: β > 0
  "H4_vs_H2" = "$H_1: \\beta$ $\\neq$ 0 vs. $H_0: \\beta$ = 0"
)

custom_y_labels <- c(
  "sidea" = "MID Side A",
  "sideb" = "MID Side B",
  "polyarchy" = "Polyarchy",
  "gdp" = "GDP cap (log)",
  "pop" = "Population (log)",
  "common_law" = "Common Law",
  "milex" = "Milex cap"
)


global_min <- min(c(m1.bf.new_long$log_BF, m2.bf.new_long$log_BF,
                    m3.bf.new_long$log_BF, m4.bf.new_long$log_BF,
                    m5.bf.new_long$log_BF, m6.bf.new_long$log_BF), na.rm = TRUE)
global_max <- max(c(m1.bf.new_long$log_BF, m2.bf.new_long$log_BF,
                    m3.bf.new_long$log_BF, m4.bf.new_long$log_BF,
                    m5.bf.new_long$log_BF, m6.bf.new_long$log_BF), na.rm = TRUE)

global_min <- global_min-.1 ## to fix the label in the figure


# Define shared fill scale
shared_fill_scale <- scale_fill_gradient2(
  low = "blue",
  mid = "white",
  high = "red",
  midpoint = 0,
  name = "log10(BF)",
  limits = c(global_min, global_max)
)

m1.new.plot <- ggplot(m1.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "",##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Overall, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m2.new.plot <- ggplot(m2.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
    geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Overall, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m3.new.plot <- ggplot(m3.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Force, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m4.new.plot <- ggplot(m4.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "Force, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m5.new.plot <- ggplot(m5.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "No force, from 1991", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )


m6.new.plot <- ggplot(m6.bf.new_long, aes(x = Comparison, y = Variable.fct, fill = log_BF)) +
  geom_tile(color = "white", linewidth = 0.5) +
    shared_fill_scale +  # Use shared scale
  geom_text(aes(label = round(log_BF, 2)), color = "black", size = 4, fontface="bold") +
  # Custom x-axis labels here:
    scale_x_discrete(name = "", ##"Hypothesis Comparison",
                    labels = TeX(custom_labels2)) +
       # labels = parse(text = names(custom_labels))) +
    scale_y_discrete(
        labels = custom_y_labels  # Apply the label mapping
    ) +
    labs(
    title = "No force, from 1945", ##"Pairwise Bayes Factor Comparisons",
    y = "" ##"Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size=14, angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5, size = 18),
    axis.text.y = element_text(size = 14)
  )

## knitr BF_2TAIL_FIGURE
BF.plot <- grid.arrange(m1.new.plot, m3.new.plot, m5.new.plot,
             m2.new.plot, m4.new.plot, m6.new.plot, 
             nrow = 2)

ggsave("gc-jg-project/kampala/paper/figures/figure3.png", plot = BF.plot, width = 16, height = 10, device = "png")


bf.1.sidea.1991 <- round(m1.bf_long[1,4],1)
bf.2.sidea.1991 <- round(m1.bf_long[2,4],1)
bf.3.sidea.1991 <- round(m1.bf_long[4,4],1)

bf.1.sideb.1991 <- round(m1.bf_long[7,4],1)
bf.2.sideb.1991 <- round(m1.bf_long[8,4],1)
bf.3.sideb.1991 <- round(m1.bf_long[10,4],1)


bf.new.1.sidea.1991 <- round(m1.bf.new_long[1,4],1)
bf.new.2.sidea.1991 <- round(m1.bf.new_long[2,4],1)
bf.new.3.sidea.1991 <- round(m1.bf.new_long[15,4],1)

bf.new.1.sideb.1991 <- round(m1.bf.new_long[3,4],1)
bf.new.2.sideb.1991 <- round(m1.bf.new_long[4,4],1)
bf.new.3.sideb.1991 <- round(m1.bf.new_long[16,4],1)

bf.new.1.sidea.no.force.1991 <- round(m5.bf.new_long[1,4],1)
bf.new.2.sidea.no.force.1991 <- round(m5.bf.new_long[2,4],1)
bf.new.3.sidea.no.force.1991 <- round(m5.bf.new_long[15,4],1)



